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ABSTRACT 

In  this  paper,  we  examine  some  splitting  techniques  for  low  Mach  number 
Euler  flows.  We  point  out  shortcomings  of  some  of  the  proposed  methods  and 
suggest  an  explanation  for  their  inadequacy.  We  then  present  a  symmetric 
splitting  for  both  the  Euier  and  Navier-Stokes  equations  which  removes  the 
stiffness  of  these  equations  when  the  Mach  number  is  small.  The  splitting  is 
shown  to  be  stable. 
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INTRODUCTION 


For  many  computational  problems  in  low  speed  fluid-dynamics,  it  has  been 
customary  to  use  the  incompressible  Euler  or  Navier-Stokes  equations.  There 
are  essentially  two  reasons  for  doing  this:  there  is  one  less  variable,  since 
the  density  remains  constant,  and  the  stability  limit  is  independent  of  the 
sound  speed.  Recently,  however,  there  has  been  increased  interest  in  studying 
compressibility  effects  even  for  low  Mach  number  fluid  flows.  The  compress¬ 
ible  equations,  unfortunately,  have  stiff  coefficients  due  to  the  disparity  in 
the  magnitude  of  the  flow  velocity  and  the  speed  of  sound.  To  overcome  this 
difficulty  various  splitting  methods  have  been  proposed  to  remove  the  stiff¬ 
ness  from  the  matrix  coefficients  of  the  equations,  [3,  A,  7).  Some  of  these 
methods,  however,  have  not  performed  as  anticipated;  in  fact,  often,  for  the 
the  stipulated  stability  limits  on  the  time  step,  the  calculations  diverged. 

In  this  paper,  we  first  propose  an  explanation  for  this  behavior.  We 
give  examples  in  the  first  three  sections  which  show  that  splittings  resulting 
in  matrices  which  are  not  simultaneously  symmetrizable  (such  as  in  [7])  may  be 
ill-posed  at  the  p.d.e.  level.  Similar  results  are  presented  for  some 
explicit  numerical  schemes,  both  finite  difference  and  spectral.  Thus,  the 
intent  of  these  sections  is  to  caution  against  unrestrained  use  of  splitting 
methods. 

In  Section  IV,  we  present  a  transformation  of  variables  which  symmetrizes 
the  Euler  equations.  Under  the  assumption  of  low  Mach  number  flow,  we  are 
able  to  propose  an  efficient  splitting  technique  for  the  compressible  equa¬ 
tions.  The  resulting  algorithm,  given  both  for  the  Euler  and  Navier-Stokes' 
equations,  is  unstiff  for  the  nonlinear  field,  and  the  other  split  operators 
are  linear  and  may  therefore  be  solved  implicitly  with  ease.  (The  implicit- 
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ness  is  necessary  to  overcome  the  stiffness  which  was  transferred  into  the 
linear  part.)  The  total  scheme  may  be  shown  to  be  stable  under  the  less  re¬ 
strictive  time  step  of  the  nonlinear  part.  In  a  future  paper,  we  intend  to 
present  computational  results  for  our  proposed  algorithm. 


1.  A  MODEL  PROBLEM 

Consider  the  initial  value  problem  for  the  following  symmetric  hyperbolic 


system 


\  ■  0  •  <J  !)  C) 


*  A  w 


(1.1) 


0  is  a  real  number,  | B 1  >  1.  The  eigenvalues  of  A  are 


u  j ( A)  =  1+0,  U2(A)  =  1  “  B, 


and  therefore  an  explicit  scheme  will  have  the  CFL  condition 


At  < 


const 


Ax. 


(1.2) 


For  example,  the  Lax-Wendroff  scheme 


n+1 


+  !£*  (wn 
2Ax  V  1+1 


n  \ 

- 


A2(At)2 
2  (Ax) 


(w 


.1  +  1 


-  2w, 


Wi-1) 


(1.3) 


Is  stable  under  the  condition 


£( 1  +  ' 
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Suppose  now  that  one  attempts  to  advance  the  solution  of  (1.1),  equation  by 
equation,  rather  than  to  use  the  form  of  the  system  as  in  (1.3).  This  amounts 
to  splitting  the  matrix  A  into  the  sum  of  two  matrices  B  and  C 


/I 

b\ 

+  1° 

0) 

'0 

0) 

+  u 

1/ 

(1.4) 


and  advancing  the  solution  by  using  first  the  equation 

(1.5) 

(1.6) 

where  the  initial  value  of  (1.6)  at  every  time  step  is  the  value  of  w^1^ 
obtained  after  advancing  (1.5)  one  time  step.  This  procedure  yields  a  scheme 
which  is  first  order  in  time  and  second  order  in  space.  We  note  that  the  sys¬ 
tems  defined  in  (1.5)  and  (1.6)  are  strictly  hyperbolic  and  hence  well- 
posed.  The  eigenvalues  of  B  and  C  are  0  and  1,  and  therefore  the  Lax- 
Wendroff  scheme  for  (1.5)  and  (1.6)  separately  will  be  stable  under  the  condi¬ 
tion 

£1  <  1  (1.7) 


and  then 


w 


(1) 


B  w 


(1) 


.(2) 


C  w 


(2) 


allowing  a  time  step  much  larger  than  the  one  allowed  by  (1.2)  if  0  is  a 
large  number.  However,  even  if  a  numerical  method  is  stable  for  (1.5)  and 

(1.6)  separately,  it  need  not  be  stable  for  the  combination  of  (1.5)  and 

(1.6) .  In  fact,  consider  the  Lax-Wendroff  scheme  for  (1.5)  and  (1.6).  The 

amplif ication  matrix  G  of  the  combined  scheme  is  given  by 
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/g(0  +  e2(l  -  g(c))2  -6g(c)(l  -  g(O) 

\  -8(1  -  g(O)  g(C) 


(1.7) 


where 


5 


=  sin 


kAx 

2 


g(S)  =  1  -  2xV  +  2iXc/l  -  t2 


X 


At 

Ax 


We  will  show  that  the  eigenvalues  of  G  are  greater  than  one  in  modulus  for 
any  X,  and  thus  the  combined  scheme  is  unconditionally  unstable.  To  do 
that  we  Look  at  the  mode  £  =  1 : 


GU  =  1) 


■( 


2  2  4 

1  -  2X  +  48  X 


-2BX 


-8(1  -  2X2)2X2 


1  -  2X 


The  eigenvalues  of  G  p*  are  given  by 


=  1  -  2X  +  26  X 


2  -  ~2-4  *  JF* 


8~X  +  1  -  2X2  28 X 2 


The  scheme  is  clearly  unstable  for 


since  In  this  case  y  +  >  1  for  8  >  1  and  u_  >  1  for  8  <  1.  It 
is  also  easily  verified  that  p  +  >  1  for  6  >  I  for  any  X.  Thus,  the 
splitting  (1.5)  -  (1.6)  is  the  wrong  way  of  splitting. 

Perhaps  a  deeper  Insight  is  obtained  if  we  Fourier  transform  (1.1), 
(1.5),  and  (1.6).  The  solution  operator  for  (1.1)  in  Fourier  space  is 
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E(w  ,At) 


iAwt 

e 


where  u  is  the  dual  Fourier  variable. 

The  solution  operator  for  the  split  scheme  (1.5),  (1.6)  over  one  time 
step  is 

iBwAt  iCwAt 

S(w  ,At)  =  e  e  .  (1 .8) 


For  every  fixed  w 

S(ai,At)t/AC  -►  E(t). 

However,  since  C2  =  C,  B2  =  B,  an  expansion  of  the  right  ^->ncj  side  of  (l.R) 
shows  that 

S(w ,A  t )  =  fl  +  B(eiwAt  -  1))|I  +  C(elw&t  -  1)]. 

If  we  put  Atm  =  it,  we  get 

/  4B2  -  1 

S(io  ,At)  =  (I  -  2B)  ( I  -  2C)  =  ( 

\  -2B 

and  for  any  1 6 1  >  1  S(o>,At)  has  eigenvalues  larger 
trates  the  instability. 


than  1.  This  illus- 


2.  THE  ISENTROPIC  EULER  EQUATIONS 

The  isentropic  Euler  equations  in  one  space  dimension  may  be  written  as 


w 


t 


A  w 


(2.1) 
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where  u  is  the  velocity,  p  is  the  pressure,  p  is  the  density,  and  y 
is  the  adiabatic  constant  of  the  fluid.  The  normalized  equation  of  state  for 
the  fluid  is 

p  =  PY.  (2.2) 


The  eigenvalues  of  the  matrix  A  in  (2.1)  are  u  -  c  and  u  +  c,  where 
c  =  is  the  sound  speed.  Thus,  if  we  were  to  solve  (2.1)  by  an  explicit 

difference  scheme,  we  would  have  to  impose  a  CFL  condition  of  the  form 


At  (  const 
Ax  —  |uj+c 


(2.3) 


We  wish  to  study  (2.1)  in  the  low  mach  number  regime  so  that  p  =  p^, 
where  is  the  base  flow  density.  We  define 


P 


Then  | € |  <<  1.  Using  (2.2)  we  conclude 


P 


-  Y€ Pg  +  0(62), 


where  pq  is  the  base  flow  pressure. 

One  possible  splitting  for  (2.1)  [7],  is  to  write  A  as  the  sum  of  two 
matrices  Aj  and  A2  as  follows 


"0 

■'p0 

U 

1/p  -  1/p0 

-Ypo 

0 

-Y(P-Pq) 

U 

(2.4) 
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We  then  advance  the  solution  of  (2.1)  by  first  using  the  equation 


t  lx’ 


(2.5) 


and  then  the  equation 


(2)  A  (2) 

wt  =  A2Wx  ' 


(2.6) 


Since  Aj  is  a  constant  matrix,  we  could  solve  (2.5)  analytically  thus  doing 
away  with  any  CFL  restriction.  The  eigenvalues  of  the  matrix  A2,  however, 
are  -  u  ±  i/y  cQ€  +  0(€  ) .  Thus,  the  splitting  (2.5  -  2.6)  is  not  a 
hyperbolic  splitting. 

To  examine  the  stability  of  the  split  scheme,  we  examine  the  Fourier 
transform  of  the  solution  operator,  S(At),  over  one  time  step.  The  Fourier 
transform  of  3  is  S: 


„  iA„u)At  iA.wAt 

S(w,At)  =  e  e 


(2.7) 


Let  a  =  Cq  uAt,  and  6  =  /y  cq  €  wAt. 


After  some  computation,  we  obtain 


iAjWA  t 


-isina/ CqPq 


“icoposlna 


To  first  order  in  € ,  we  may  write  A2  as 


-r  p0€ 
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We  then  have 


lA^oiAt 

e 


iuwAt 

e 


coshB 

-iA  CgppSinhB 


-isinhB/A  cgP0‘ 
coshB 


Hence , 


rcoshBcosa  - 


S(ui,At)  =  e 


iuwAt 


sinhBsina 

/7 


-i 


copo 


(coshBsina  +  sinhBcosa)-) 


Lip0cn(/Y  sinhBcosa  -  coshBsina)  A  sinhBsina  +  coshBcosa 


The  eigenvalues  of  S(io,At)  are  roots  of  the  polynomial 

p(X)  =  X^  -  (2cosh8cosa  +  — — —  sinhBsina)X  +  1=0. 

A 

By  Miller's  criterion  [8]  the  roots  of  p(X)  are  inside  the  unit  disc  if 
and  only  if 

0  =  ] coshBcosa  +  -  sinhBsina |  <  1. 

2A 

If  we  let  a  =  Cq  coAt  =  tt  ,  then  8  >  I,  whenever  0  >  0.  Hence,  at  least 
one  of  the  eigenvalues  of  S(ui,At)  lies  outsidt  the  unit  disc  in  a  neigh¬ 
borhood  of  a  =  n  . 

If  we  were  to  solve  (2.5)  -  (2.6)  using  a  pseudospect ra 1  difference 
scheme,  we  would  have  to  impose  a  CFL  restriction  of 
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difference  scheme  for  (3.1)  would  have  a  CFL  restriction  of  the  form 


At  .  const 

a  __  ' 


Ax  —  TuTTc 


(3.2) 


One  possible  splitting  for  (3.1)  could  be  obtained  by  writing  A  as 


'0 

1 

o' 

0 

0 

0  ' 

A  =  - 

0 

0 

1 

2 

-u 

2u 

0 

.0 

2 

co 

0. 

2 

-c  u 

2  2 
c  -c0 

u_ 

where  c0  is  the  sound  speed  of  the  base  flow, 
first  the  equation 


Aj  w 


(1) 

x 


and  then 


(2)  .  (2) 
w  '  =  Aw  , 
t  2  x 


We  then  solve  (3.1)  by  using 

(3.4) 

(3.5) 


The  advantage  of  such  a  splitting,  it  would  seem,  is  that  since  Aj  is  a 
constant  matrix  we  can  obtain  an  analytical  solution  of  (3.4)  without  any 
restriction  on  the  time  step.  Further,  since  the  eigenvalues  of  are  0, 

-u,  and  -2u,  we  can  solve  (3.5)  by  a  difference  scheme  with  a  large  CFL  con¬ 
dition  of  the  form 


At  .  const 
Ax  —  2 |u| 


(3.6) 


We  examine  the  Fourier  transform  of  the  solution  operator  S(At)  over  one 
„  iA^wAt  iA^uiAt 

time  step.  Then  S(w,At)  =  e  e 


Let 


a 
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n  = 


2  2 
e  -  c0 


We  choose  u  *  0  and  n  >0.  Then 


iAjtoAt 


0 

Lo 


islna 


cosa 


-Ic^slna 


1+cosa  — i 
2 

c0 

Islna 


0 

cosa 


and 


lA^uAt 


1 

0 

Li 


o 

l 

-inc0a 


0 
0 
1 J 


Hence 


S(ao  ,At) 


islna 

C0 


cosa 


-Knc^a  cosa  +  c^slna) 


1+cosa 
2 

C0 

Islna 
C0 

(-naslna  +  cosa)J 


The  eigenvalues  of  S(w,At)  are  1,  and  the  roots  of  the  polynomial 


p  (A)  =  A  -(2cosa  -  naslna)  .+  1  =  0. 


By  Miller's  criterion  the  roots  of  p( A )  are  inside  the  unit  disc  if  and 


only  if 


-12- 


Let  a  =  — it  +  <5  . 

'  Then  e  =  1  +  nw6  +  0(6 2) . 

1  Hence  S(ai,At)  has  at  least  one  root  outside  the  unit  disc  near  a  = 

CpWAt  =  it.  Thus,  this  proposed  splitting,  once  more,  has  undesirable  proper¬ 
ties.  If  we  were  to  solve  (3.4)  -  (3.5)  using  a  pseudospectral  difference 

•  scheme,  we  would  have  to  impose  a  very  restrictive  CFL  condition  of  the  form 

At  „  const 
Ax-  cQ 

Gustafsson  and  Guerra  [3]  showed  how  to  split  (2.1)  in  a  way  that  avoids 

^  the  pitfalls  pointed  out  above.  The  main  idea  in  their  work  was  to  obtain  two 

symmetric  split  operators.  This,  of  course,  is  harder  to  do  for  more  compli¬ 
cated  systems.  In  the  following  sections,  we  generalize  this  approach  to  the 
problem  of  obtaining  split  operators  which  are  simultaneously  symmetr izable  in 
the  case  of  the  full  Euler  and  Navier-Stokes  equations. 


4.  CORRECT  SPLITTING  FOR  THE  EULER  EQUATIONS 

In  the  preceding  sections,  we  gave  examples  of  "natural"  splitting  proce¬ 
dures  which  led  either  to  Instabilities  or  to  stability  conditions  which  at 
best  did  not  represent  an  improvement  over  the  original  ones.  A  common  fea¬ 
ture  of  those  split  operators  was  that  they  were  not  simultaneously  symme- 
trizabie. 

To  avoid  the  dangers  pointed  out  by  these  examples,  we  propose  to  remove 
the  stiffness  of  a  given  stable  symmetric  operator  by  instituting  a  splitting 
procedure  such  that  all  the  split-off  operators  are  simultaneously  symme- 
trizable.  If  each  of  these  new  operators  is  discretized  in  a  stable  manner, 
then  the  overall  scheme  will  remain  stable. 
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A  prescription  for  a  general  operator  achieving  this  goal  Is  not  known  to 
us.  We  would  like,  however,  to  suggest  such  a  procedure  for  compressible,  low 
Mach  number  flows  governed  by  either  the  Euler  or  the  Navier-Stokes  equa¬ 
tions.  These  systems  are  chosen  in  view  of  the  "counter-examples"  given  In 
Section  3.  The  Euler  equations  may  be  symmetrized  nonlinearly  by  using 
"entropy-variables"  [5,  6].  The  system  thus  obtained  is  of  the  form 


pq. 


3 

1 

1=1 


a.  q 

1  X  . 


=  o 


(4.1) 


where  P  and  the  A^'s  are  symmetric  maLrix  functions  of  the  vector  q.  The 
premultiplying  matrix  P  is  usually  non-sparse,  and  hence  It  is  not  clear  how 
to  remove  the  stiffness  (If  there  is  any)  from  the  A^'s.  In  the  Euler  equa¬ 
tions,  it  is  well  known  that  the  eigenvalues  of  Aj  are  u,  u  +  c,  u  -  c 
where  u  is  the  x-componenl  of  velocity  and  c  Is  the  speed  of  sound.  At 
low  Mach  number  flows,  u  <<  c  everywhere;  hence,  a  Von-Neumann  like  stability 
condition 


At  < 


Ax 

M+c 


(4.2) 


gives  an  over-restricted  condition.  In  this  sense,  the  system  is  stiff  (see 
Sections  2  and  3). 

Our  approach  is  motivated  by  previous  results  [1]  valid  for  the 
linearized  frozen  coefficient  case. 

Consider  the  Euler  equations  for  a  gas  In  their  nonconservative  form  in 
two-space  dimensions  (the  three-dimensional  case  follows  directly  from  the 


results  of  this  section): 
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>|<> 

+ 

b^  = 

0 

(4.3) 

at 

3x  3y 

where 

V 

is 

the  column 

vector  (p ,u , v, p) 

and  the  coefficient 

matrices 

are  given 

by 

u 

P 

0 

0 

V 

0  p 

0 

A 

= 

0 

u 

0 

1/p 

,  B  = 

0 

v  0 

0 

,  (4.4) 

0 

0 

u 

0 

0 

0  v 

1/p 

.0 

YP 

0 

u  . 

.0 

0  YP 

0  . 

where 

P  1 

,u,v,p 

and 

Y 

are, 

respectively , 

the 

density,  the 

velocity  com- 

ponents 

in  the 

A 

X 

and 

y 

directions , 

the 

pressure  and 

the 

ratio  of 

specific  heats  at  constant  pressure  and  volume 
nondimenslonalize  the  quantities  in  (4.3)  as  follows: 


(,  .  Cp/Cv>. 


Next , 


tU 

“>  x  y  p 

1  =  ~T~  ,x=T,y=T’p=p~ 


u  v 

u  =  —  ,  v  =  —  ,  p 
u  v 


(4.6) 


P  u 

00  oo 


where  the  subscript  «  indicates  free  stream  condition  and  i.  is  a 
reference  length.  Equations  (4.3)  and  (4.4)  then  retain  the  same  form  exactly 
with  the  superscript  *  removed.  In  particular,  the  dimensionless  speed  of 
sound  retains  its  functional  form,  i.e., 


c  =  —  =  /yp/p  • 


(4.7) 


We  now  propose  the  following  change  of  variables: 
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(4. ft) 


where  Cj  is  a  constant  to  be  specified  later.  One  may  then  cast  (4.3)  in 
the  form  of  (4.1)  where: 


(4.10) 


0 


(4.11) 
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are  symmetric  hyperbolic. 

We  now  wish  to  motivate  the  way  in  which  the  operators  in  (4.11)  (Aj  ^ 

an(j  A„  l—  )  are  split.  The  starting  point  is  the  fact  that  we  are  here 
2  3y 

interested  in  low  Mach  number  flow.  Such  flows  are  characterized  by  two 

2  2  2 

facts:  the  first  is  that  c  »  u  +  v  everywhere;  secondly,  for  reason¬ 


able  initial  conditions 


2/  \  2 
c  (x,v,t)  -  c 


«  1. 


(4.12) 


For  example,  at  steady  state 


c^(x,y)  -  c“  T  -  T 
* J  °°  ^  st  « 


(4.13) 


where  T  ^  is  the  stagnation  temperature,  M^  is  the  free  stream  Mach 

number;  hence,  for  low  Mach  numbers  (4.12)  is  valid.  In  view  of  the  above,  we 

choose  c,  =  c  ,  and  we  rewrite  (4.11)  as  follows 
1  °° 


PqL  +  (Rj  +  Sj)qx  +  (R2  +  S2)qy  =  0 


(4.14) 


where 


The  four  eigenvalues  of  P”*Rj  are 


X  =  u,  u,  u  ±  (c  -  c^Ml  + 


It  Is  clear  from  (A. 11)  that 
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c  -  c  <  ( y  —  1  )M 

CO  _  QO 

while 

—  -  0(1). 

c 

Thus,  none  of  the  eigenvalues  gets  to  be  large  unlike  the  original  unspill 
scheme  which  had  eigenvalues  u,  u,  u  ±  c  (recall  that  In  our  case  u  “  0(1) 
while  c  “  1/M  .) 

OO 

Next  we  notice  that  Sj  and  S2  are  constant  matrices.  A  difficulty 
remains  however  In  the  nonlinear  element  of  P.  We  shall  deal  with  this  as 
the  method  of  solution  Is  presented. 

Step  I  in  the  solution  of  (4.12)  is  to  numerically  advance 

Pqt  +  Vx  +  R2qy  =  0  (4*15) 

by  one  time  step.  We  have  just  demonstrated  that  stability  criterion  for 

(4.15)  is  not  stiffly  restricted.  in  fact,  for  most  explicit  schemes,  to 
within  a  constant  of  order  unitv, 

^  ,  ;  Ax  Ay  1 

-  x ,  y  '  Ju  |  + 1  c-c^  |"  ’  |v|  +  |c-cj  1  ’ 

as  compared  with  (4.2).  The  gain  is  obvious.  Step  II  in  the  procedure  Is  to 
solve 

Pqt  +  Slqx  +  S2qy  =  °*  (4*lb) 

The  Initial  conditions  for  (4.16)  are  given  by  the  solution  to  (4.15)  at 

t  =  At.  Notice  that  while  Sj  and  S2  are  constant  matrices,  P  has  the 
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2  2 

nonlinear  element  c  / c  .  This  means  that  removal  of  the  stricter  time-step 
due  to  Sj  and  S2  (At  £  const.  Ax/c^)  cannot  be  done  easily  via  implicit 
method  implementation  of  (4.16).  To  overcome  this  difficulty,  we  split  Sj 
and  S2  as  follows: 


where 


and 


S  -  sT  +  s11 

bl  ■  bl  bl 


S  =  sr  +  sTI 
b2  S2  b2 


'0 

c 

00 

0 

rO 

0 

0 

0 

— 

o' 

c 

00 

0 

0 

cii 

0 

0 

0 

,Y-1 

/ -  C 

0 

Y  00 

✓7 

si 

0 

0 

0 

0 

0 

.  0 

0 

0 

0 

0 

O  O 

l_ 

-0 

,r-l 
/ -  c 

y  00 

0 

0 

(4.17) 
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0 

✓7 

0 

0 

0 

0 

0 

0 

0 

0 

CII 

S2  = 

c 

OO 

0 

0 

L 

0 

0 

0 

/EE  c 

— 

0 

Y  CO 

V  — 

-  c  0 

_0 

0 

0 

0. 

0 

0 

/— 

Y 

00  J 

(4.18) 


Thus,  we  replace  (4.16)  with  the  sequence 


Pq,  +  S*q  +  S*q  =0 
t  lx  2  y 


(4.19) 


Pq  +  +  S^q 

Mt  1  Mx  2  Hv 


0 


(4.20) 
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Note  that  in  (4.19),  ct  is  identically  zero,  and  so  over  that  time  step  we 
take  c  =  c(x,y)  from  (4.15).  The  rest  of  (4.19)  is  therefore  linear  (be- 
cause  =  c  (x,y)  is  known)  and  can  be  solved  implicitly  with  relative 

ease.  Alternatively,  the  first  three  equations  in  (4.19)  may  be  combined  into 
a  variable  coefficient  wave  equation  for  in  p,  namely: 

3 2  c«  ? 

— r-  (in  p)  =  — r -  V  (in  p);  (4.21) 

3t  c  (x,y) 

u  and  v  are  then  obtained  directly  from  the  middle  two  equations  of 
(4.19).  In  (4.20)  it  Is  in  p  that  does  not  change  over  the  time  step. 
The  rest  of  the  system  is  linear  with  constant  coefficients  and  may  also  be 
cast  Into  a  wave  equation  for  c: 


Y  — 1 
Y 


2 

c 

00 


V2c , 


(4.22) 


and  again  u  and  v  are  found  directly  from  the  middle  two  equations  of 

(4.20). 

This  completes  the  splitting  method  for  the  Euler  equation.  The  temporal 
and  spatial  discretization  depends  on  the  particular  problem.  Straight 
splitting  as  described  here  will  result  in  only  first  order  accuracy  in 
time.  Alternating  the  order  of  solving  between  (4.15)  +  (4.19)  +  (4.20) 
to  (4.20)  +  (4.19)  ♦  (4.15)  will  yield  second  order  in  time,  see  [2). 
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5.  EXTENSION  TO  THE  NAVIER-STOKES  EQUATIONS  (N.  S.  CASE) 

The  Navler-Stokes  equations  may  be  written  as: 


Pq  +  A  q  +  A  q  =  B  q  +  B  q  +  Cq 
t  l  x  2  y  1  xx  2  yy  ^ 


+  F,  + 
xy  1  2 


where  q,  P,  Aj,  and  A2  are  as  defined  in  the  preceding  section, 
quantities  on  the  right  hand  side  are  given  by: 


B ,  =  ~ 


1_ 

Re 


4_ 

3 

0 

0 


0 

0 

1 

0 


0 

0 

0 
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2Prp 


=  rr- 
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Re 


0 

0 

4_ 
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2Prp 


C  =  Re 


0 

0 

J_ 

3 

0 


0 

1 

T 

0 

0 


F1  =  — 
1  /T-T 


0 

0 

0 

1 

P  R 
'  r  e 


i  <^> 

P  c 


(5.3) 


F2  " 


1 


7VT7-1 ) 


JL 

Re 


0 

0 

0 

pc 


where  the  dimensionless  viscous  production  function  $  is  given  by 


~  +  v  +  2[u^  +  v^ ]  +  [u  +  v  ]^. 

J  X  y  X  y  V  X 


(5.1) 

The 

(5.2) 

(5.3) 

(5.4) 

(5.6) 


(5.7) 
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We  can  now  describe  the  solution  method:  after  obtaining  the  "hyperbolic" 

solution  (see  equations  (4.15)  to  (4.22)),  we  go  through  the  following  steps: 

1)  From  (5.2)  to  (f.6).  It  follows  that  =  0,  l.e.,  during  the  "viscous" 

integration  p  =  p(x,y). 

2)  multiply  equation  (5.1)  by  the  matrix 

'0  0  0  0  ' 

0  10  0 
0  0  1  0 

.0  0  0  0  . 


The  resulting  two  "viscosity  split"  equations  for  u  and  v  are 

(5.8) 

(5.9) 

These  may  be  easily  solved  implicitly  since  they  art  linear  p.d.e.'s  with  con¬ 
stant  coefficients. 


u  *  -i—  [  (ir  u  +  u  )-*-i-u  ] 

t  Re  3  xx  yy  3  xy 


Vt  Tfe  ^Vxx 


+  7  v 


+  4  v  ], 

xy  3  xy 


(3)  The  last  step  is  to  solve  the  viscous  part  of  the  energy  equation  which 
may  be  cast  In  the  form: 


3 

5T 


(c2) 


y 

SPrRep 


[V2 (c2) ]  + 


1  4> 

7e“p  * 


(5.10) 
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Note  that  $/p  Is  a  function  of  the  squares  of  ux,  uy,  vx,  vy,  and 
p(x,y)  only  and  may  therefore  be  taken  as  known  from  the  previous  step. 

Equation  (5.10)  is  then  a  scalar  linear  Inhomogeneous  heat  equation  which 
again  may  be  easily  solved  by  implicit  methods. 

Note  that  all  the  operators  in  (5.8)  and  (5.9)  may  be  taken  to  be 
stable  (i.e.,  11*11  <  1)  in  •  In  addition,  because  c^  >  0,  F2  >  0 
(5.10)  is  also  stabillzable  under  the  Lj  norm  for  c^;  this  assures  the 
L2  stability  for  q^ . 

Notice  the  total  algorithm  (4.15)  -*■  (4.22)  +  (5.8)  +  (5.10)  may  be 
run  partly  in  parallel  thus  enhancing  its  efficiency  beyond  the  removal  of  the 
stiffness.  Schematically,  the  tree  of  calculation  may  be  shown  as  follows: 


(eq .  4.21) 


nonlinear  part  of  Euler  (eq.  4.15) 


CD 

(eq.  1 4 . 2  2 ) 


(eqs.  5.8, 5.9) 


~sr" 

(eq .  5.10) 


Thus,  If  parallel  processors  are  available,  we  run  only  three  computations 
Instead  of  five. 
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